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ABSTRACT 


There are savaral important causes which can make the nuclear 
power plant unstable and ultimately lead to an accident . One of 
such an important causes (faults) is small break loss of coolant 
accident ( SBLOCA ) . Slow but steady inventory loss from the 
steam generator (SC) coolant system will degrade the performance 
in the long run . It is useful to detect the SBLOCA at proper 
time . The fault detection can be done in several ways . Out of 
which General Likelihood Ratio ( CLR ) test is the most useful one 
Hitherto , the CLR method has been appl ied by other researcher 
only to ident i f y the sensor fai lures in nuclear power plants . In 
this work we have attested to apply the CLR method to detect 
SBLOCA in the steam generator portion of the primary heat trans- 
port system ( PHT ) . We developed a simpl if ied state space vari- 
able model of steam generator and tested the utility of the CLR 
method . We also developed a simple digital control system to con- 
trol the outlet pressure of the primary coolant by manipulating 
the flow rate of the coolant . And level of water is control by 
manipulating the flow rate of water . The result shows that a leak 
of magnitude as low as 0.25X of the total flow rate of primary 
coolant can be detected 90 t imes out of 1 00 times within SO . 0 
seconds (approx) . Thus we can conclude that SBLOCA of the PHT 
loop can be detected by CLR method . Further study should be done 
by modelling the complete PHT loop and also the other prob- 


able accidents . 
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INTRODUCTION 


1.1 PROBLEM CONSIDERATION 

Among the different nuclear power reactors pressurized heavy water 
reactors ( PHURs ) have few more advantages over others because of 
it's simplicity in operation and use of natural uranium ( 0.7 X U- 
£35 ) as fuel . In India , nuclear plants are mainly based on PHWR; 
at present four PHWR are in operation . Nuclear power plants ( PHUR 
type ) have several units which are inter-connected in a close loop 
and all PHUR type plants have common primary heat transport ( PHT ) 
loop which is shown in the fig 1.1 . Our present problem is to detect 

leak ( only due to small break ) in the loop using a model based on 
statistical method , called GENERAL LIKELIHOOD RATIO ( GLR ) method 

1 . £ REVIEW OF THE PAST WORK 

The GLR method was developed by Uillsky ft Jones in 1974 . Though i 
is used to detect faults in many fields , it's applicat ion in nuclear 
field is very few . One such application is fault diagonisis in pressu 
rizer ( ref. 14) .Another attempt is made to detect SBLOCA in the heav 
water reactor ( ref. 13 ) and it is reported that a leak of magnitude 
0 . £5 X can be detected with 98X accuracy . 

1.3 THE PRESENT WORK 

PHT loop contains several units . All the units should be modelled 
to detect leak in the loop . Only steam generator is considered for 
present work . In this present work a simplified model is developed 
for normal operation as well as in the presence of a leak . It is assu? 
med that the plant is operating around a normal steady state and thus 
the model will not be valid for start-up ft shutdown period . 




gg LIGHT WATER MODERATOR HEAT EXCHANGER 



LIGHT WATER STEAM 
HEAVY WATER COOLANT 


HEAVY WATER MODERATOR 


FIG.] A SIMPLIFIED DIAGRAM OF PHT LOOP 






CHAPTER - £ 


THE CLR METHOD 


GENERAL LIKELIHOOD RATIO ( CLR ) is a statistical procedure to 
detect and identify a fault , if any in the system . It was first 
introduced by Uillsky ft Jones ( ref. 16 ) in 1974 . Let us re- 
view the method in brief . Further details are available in refer. 

( 16 , 19 , 20 , £1 ) 

S.1 PRINCIPLE OF THE 6LR METHOD 

Let us first consider the CLR method pictorially . Suppose we 
are given the pattern of differnt measurements with respect to time 
for different faults . Every fault has specific pattern . This stan- 
dard pattern can be matched with the plant data obtained and the 
fault is declared when the trend of data closely matches with the 
pattern of a particular fault . Fig 2.1 explains this clearly . The 
fig 2. la and 2.1b shows the analytically generated pattern for 
temperature and pressure for two different faults 1 ft 2 .In 
fig 2.1c we plotted schematically the values of the same variables 
obtained from the plant . This is showing that the pattern generated 
by the plant data is matching with the analytically generated pattern 
of fault- 1 . 

In the above approach two important questions that need to be 
answered are ; 

(1) HOU TO OBTAIN THE STANDARD PATTERNS FOR DIFFERENT FAULTS 

(2) HOU TO MAKE THE DICISION THAT A PARTICULAR FAULT PATTERN 
BEST FITS THE MEASURED DATA . 

The CLR method answers to both the questions , which is stated as 
follows ; 

In order to obtain the standard patterns for different faults 
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th* GLR method needs a linear state space model of the process when 
the fault is present , that is , a model for every possible faults 
that needs to be detected in the plant . Based on this model , the 
GLR method analytically computes the expected pattern of measurement 
for the fault in a recursive fashion .The mathematics of this recur- 
sive computation is described later . These standard patterns for 
each fault , thus obtained , are stored for subsequent comparision 
with measured data . 

For obtaining the goodness of fit between a standard fault 
pattern and measured data the GLR method uses a statistical test 
called LIKELIHOOD RATIO TEST . A statistical test is necessary 
since we assume both the process model and the measurements are 
corrupted by random noise due to f luctuat ion in some of the 
inputs . The 1 i kel ihood ratio test can be viewed as a generalization 
of the standard least square fitting procedure in which the measure 
ments are weighted inversely according to their variances . The 
likelihood ratio for any fault is the ratio of two residuals , of 
the weighted least square fit , the residuals obtained when the 
measured data is fitted with the pattern for no fault to the residu- 
als obtained when the same measured data is fitted with the standard 
pattern for the fault . Under the assumt ion that all random noise 
are normaly distributed the 1 i kel ihood ratio is also equivalent to 
the rat io of two probabilities , the probabi 1 ity that a fault has 
occurred to the probability that no fault occured for the given 
measurements . Thus , if this ratio exceeds a pre specified limit , 
we can concluded that the fault has occurred . Among different faults 
, the fault corresponding to the maximum likelihood ratio will be 
identified as a fault if it's maximum likelihood ratio exceeds the 
pre-specif ied limit . 

Although we described the procedure in term of measured data , 
the actual comparision is made using what are called " INNOVATION " 
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or the difference between Measured and predicted values which is 
described below . 


22 BRIEF MATHEMATICAL REVIEW OF THE METHOD 

Let us consider a general state space representation of the 
Modelled equations as follows : 


X 

- k + 1 


A 


X + B . U + D . U 

- k - - k - - k 


U « C . Y 

- k - - k 


Y - H . X + V 

- k - - k - k 


where the vectors : 


X is deviation of state variables at time k 

- k 


U is deviation of the Manipulated variables 

- k 


U is deviation of noise at tiae k 

- k 


V is of the noise in aeasureaents at tiae k 

- k 


2.1 


2.2 


2.3 



Y is deviation of controlled variables at tine k 

- k 

and matrices : 

A is the state matrix 

C is the control matrix 

H is measurement matrix 

D is the noise matrix 

Now the GLR method can be used in in the following sequence ; 

- Kalman filter is used to estimate variables ( X ) 

- k + 1 

Pattern of innovat ions are computed 

- GLR test is performed to compare with pre-determined 
patterns 

If fault is present , it's magnitude is declared . 

Mathematical description of various steps are as follows : 

2 . S-A THE KALMAN FILTER 

The optimal estimate of the state variables and the error covariance 
matrix are obtained by the use of this filter . The Kalmann filter is 
basically a least square estimator , and it is minimizing error between 
the actual plant state vector and the state estimate . Mathematics 
behind the generation of estimate by Kalmann filter can be described 


as follows : 
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Lot X is conditional state estimate of the state 

- k/k-1 

variables for given measurements upto and including time ( k-1 ) . 
Corresponding to these est imates the co-varience matrix of the error 
in these estimate is defined as : 


P « cov C ( X 

- k/k-1 - k/k-1 


X ) / Y 3 2.4 

- k - k-1 


P « cov l ( X - X )/ Y 3 2.5 

- k/k - k/k - k - k 


The assumtion of the estimates are given in reference ( 21 ) . The 

computation of the estimates are as follows ; 


X -X + K CY - H . X 3... 2.6 

- k/k - k/k-1 - k - k - - k/k-1 


where , 


X « A . X +B.U 2.7 

- k/k-1 - - k-1 /k-1 - - k-1 


U = C . X 2.8 

-k-1 - - k-1 /k-1 


K is Kalman gain matrix 

- k 
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/ -1 

K “ P HEH.P . H + 03 2.9 

- k - k/k-1 - - k/k-1 - k 


/ / 

P ■ A . P .A + D.R.D 2.10 

- k/k-1 - - k-1/k-l - - 


P « C I - K . H 3 P 2.1 1 

~ k/k - - k - - k/k-1 


/V 

Initial conditions : X « E C X 3 2.12 

- 0/0 - 0 


Now innovations are given by 


a - y H . X 2.13 

- k - k - - k/k-1 

It can be used to identify gross error . 

Here E - denotes expectation values . 

cov. - denotes cross correlation between the 
parameters . 

/ (prime) - denotes transpose of the given matrix 

It should be mentioned that the measurements are not used for fault 

diagnosis directly because we do not know the expectation value of 

the measurement variables , y 

- k 



2.2-B 


THE INNOVATION PATTERN FOR THE FAULTS 


The expect ed pattern of innovations that will be observed after 
a particular fault occurs at tine t , can be computed analytically 
by making use of the fault aodel and Kalaann filter equations . 

But here we shall simulate the fault (only leak ) that would occu- 
re at the begining of the simulation .The fault model for leak is 
given below . 

£.2- B1 FAULT MODEL ( PROCESS LEAK ONLY ) 

A fault aodel describes the effect of the fault on the process . 
Here only one type fault ( ie, leak ) is considered . The models for 
other types of faults are given in ref. (13,16,18,19). 

PROCESS LEAK MODEL : A leak of magnitude b in the process 
unit j that occurs at time t*0 can be modelled as step chage 
in the state eqations , ie 

X « A . X +B.U + b . m +D.U 2.14 

- k + 1 - k - k ““j - k 

2.2-B2 EFFECT OF THE FAULT ON THE INNOVATION 

The state variable , at any time t*k , depends on the values 
of the variables at previous times , therefore the net effect 
of fault on the state variables and the innovations change from 
time to time . Ue can compute this effect by making use of appro 
piate fault model and the filter equations . 

The expected value of the innovations for any fault i of 
magnitude b that can occurs at time t » 0 which is less than 
the time of interest k can be obtained as ; 



E [ < ] m b C 

- k/k-1 - k , i 


wh* r * G is on* of th* signature Matrices and they 

- k, i 


have a general expression which is given as follows ; 


G -HIT A . J 3 

- k - - k - - k 


Where , 


T * A . T 

- k - k- 1 


J « A . J 

- k - k - 1 


£.17 


£.18 


Detailed derivation of these expressions are given in ref. C 13 , 
16, IS, 19). 

Thus the hypothesis of the fault ( leak ) detection can be 
formulated as : 

FOR NO LEAK 



H 

: 

E t 

a 

/ 

y 

3 = 

0 

. . * * 

£.19 


- 0 



- k 


- k-1 





FOR 

LEAK 











H 

. 

E C 

a 

/ 

y 

3 « 

b . 

G . m .... 

£ . £0 


- 1 



- k 


- k-1 



- k - 


B3 


IDENTIFICATION 

AND 

ESTIMATION 

OF 

LEAK 



With th* hypothesis given by the equations ( £.19 ) & ( £.£0 ) 
we will performed the GENERAL LIKELIHOOD RATIO ( CLR ) test to iden- 



tify leak . For this we have to calculate the likelihood ratio 
, which is defined as ; 


Pr ( X , >U ) 

- k 


£.21 


Pr ( X 


) 


Where , 


Pr ( X , 4* ) 

- k 


E exp C-.5 


N 


- k 


the probability of obtaining the 
observed innovations if a leak of 
Magnitude b is present . 


/ -1 

).Q . ( a )>] / Kpr 

-k - k - k -k 

£ . ££ 


Pr ( X 


0 ) * the probability of obtaining observed 

innovations if no leak is present . 


■ £ exp 


N 



/ 


a 

-k 


Q 

-k 


a ) ) 3 / Kpr 
-k 


Kpr 


. SNn 

( £ rt ) 


n 


Qi 


.5 

i*1 


£ . £3 


£ . £4 


4l«bC.n ....£.£5 

-k -k -i 


Q « error coverience matrix of a 
- k - k 



12 


Nn * no of sampling periods that transpire , subsequent to the 
occurrence of leak . 

Now the maximum likelihood ratio is given by ( ref. 18 ) 

T * 2 In ( X ) - Sup T 2.26 

i i 

Where , 


T - 


i 


2 

d 

i 



2.27 


d * 
i 



X 

- k 


i 



-1 

Q . G 
- k - k 


) . m 


-i 


2.28 


2.29 


The maximum value of T gives T . Using the principle of likelihood 

i 

test , the fault which is correspond to T , is identified and it's 
magnitude is estimated by the corresponding b 

i 

2.3 PROCEDURE FOR IDENTIFICATION OF LEAK USING GLR METHOD 

Leak is detected and estimated using a computer programe . 
Complete flow chart of the programe is shown in fig. 2.2 



ssDEV.IN MESD. VAR ssDEV.IN MESD. VAR 


FIG 2.1a FAULT 1 



TIME IN MINUTES 


0.200 FIG 2.1b FAULT 2 

1 
— T 

i *jl*jl* TEMP, 

j ***** PRES. 


n 



TIME IN MINUTES 



0.200 FIG. 2-1 c PLANT DATA 



0.00 10.00 ' 20.00 30.00 
TIME IN MINUTES 


FIG. 2.1 SCHEMATIC PLOTS TO 
EXPLAIN GLR METHOD 
P I CTO R! ALLY 



BEGIN 



FIG £ . 2 : FLOU CHART OF GLR METHOD ( simplified > 







CHAPTER - 3 


MATHEMATICAL MODELLING OF STEAM GENERATOR 


3.1 STEAM GENERATOR 

The steam generator that has been modelled is a vertical shell 
and U-tube evaporator. It produces a mixture of water and steam at 
saturated conditions . Centrifugal moisture separators located above 
the U-tube bundle , act so as to produces steam of required quality. 

A simplified model of the steam generator is developed with the 
objective of verifying whether SBLOCA can be detected . 

3.2 GEOMETRY OF THE MODEL 

The steam generator has two cylindrical sections of different 
cross sectional area as shown in the fig 3.1 . Lower section con- 

tains a bundle of U-tube and it has a smaller cross sectional area 
than the upper section . Pressurized heavy water, which is primary 
fluid , flows through the U-tube bunddle and water which is secon- 
dary fluid , flows through the shell. 

3 . £-A VARIABLES OF INTEREST -: 

Following variables are used in modelling the SG . 
i) STATE VARIABLES 

h - water level in shell 

DoS - density of steam produced 

AP - pressure drop in primary coolant tube 

Tol - outlet temperature of primary coolant 
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ii) MANIPULATED VARIABLES 

FI - flow rate of primary coolant 

( It is manipulated to control the level of fluid 
in primary coolant tube ) 

Fi2 - inlet flow rate of secondary fluid 

(manipulated to control the level of fluid in shell 
side of the 86) 

In addition , the quality of steam, X which is defined 
as mass ratio of dry saturated steam to the water vapour mixture 
, is used . The steam generator model is described in term of 
deviation variables , ie in term of the deviation of variables 
from their steady state values . Deviation variables are defined 

as ; 

_h « hs - h 

DoS « Do£S - Do2 

AP ■ PS - P 

Tot « TolS - Tol 

3.3 ASSUMPTION 

i) NEGLIGIBLE SENSIBLE HEAT : Secondary fluid enters the steam 
generator at it's boiling point at the corresponding pressure 
Thus the sensible heat transfer to the secondary fluid is 
zero . 

ii) SMALL DEVIATION : The process is assumed that it is operated 
around normal steady state value , ie , deviation is small . 
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iii) LINEAR TEMPERATURE DISTRIBUTION : Temper At urt variation of 
primary coolant is linear along tha langth of tha U-tube . 

iv) EQIDISTRIBUTION OF THERMAL IMBALANCE: Thermal imbalanca (B) is 
tha diffaranca batwaan tha haat content of in & output streams 
This imbalance is due to tha accumulation of haat in both tha 
primary and secondary side .It is assumed that this tharmal 
imbalanca is share equally by primary coolant and secondary 
fluid . 

3.4 DERIVATION OF THE MODEL 

3.4-A REDUCED BALANCE EQUATIONS FOR STEAM GENERATOR 

a) PRESSURE DROP EQUATION : Mechanical energy balance equation 
for primary coolant tubes 

Using usual meaning of friction factor , f , pressure drop in 

tha tubas is given by : 

2 

1 D1 . V 

AP « f . 

id £ 

Where V is tha velocity of coolant through tubas. 

If n is tha no of tubas and a - is cross sectional area then 


1 


id 



2 . 


AP - f . 


n 


3.1 



Using deviation variables and linearizing the above eqution 
about the normal steady state values we get 
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AL m f - 


1 

Id 


D1 . FIS 

pi 

n . a 


AP * C9 . FI 


3.2 


b) OVERALL MASS BALANCE FOR THE SECONDARY SIDE OF THE STEAM GENERATOR 


Rate of mass in - Rate of Mass out * Rate of Mass accuaulat ion 


d 

Fi2 . Di2 - Fo2 . Do£ ■ — C ( L - h) A Do2 + h A Di£ 1 .. 3.3 

dt 


Writing the above equation in teres of deviation variables we get 


/ 

( A Di2 


d _h d Do£ 

A Do2S ) + A (L - hS) - Die . Fi£ - Fo2S.Do£ 

dt dt 


Do2S . Fo2 is snail and is neglected 


Where teres involving product of deviation variables are beglected 


d _h_ d Do2 

al . + aE . alO Fi£ + a6 Do2 3.4 

dt dt 

c) OVERALL THERMAL ENERGY BALANCE IN THE SG 


Rate of heat in - Rate of heat out « Rate of heat accumulation in the 

primary + secondary side 


d EHP3 d CHS! 

or, thermal imbalance (B) + 

dt dt 
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U.ing th. assumption (iv) and writing th. aquation in t.r.s of d.vi 
at ion variables we get 


1 d tHP3 

2 dt 


3.6 


t d CHS] 

• 1 ■ — 

2 dt 


3.7 


Now we have to derive an expression for ACCUMUMATION OF HEAT IN PRIMARY 
AND SECONDARY SIDE 

PRIMARY SIDE : Fig 3.1 shows one U-tube of length 1 and cross sectional 
area , a . Let us consider a snail elenent of length dx at a distance x 
fron the inlet of the coolant. Now using assuntion (iii) temperature of 

the element is given by 

Tx - Til + x . (Tol - Til ) / 1 
Heat content of the small element is 
- a . dx . D1 . CpI . Tx 
Total heat content for a single tube is 

rx~l 

« ( a . D1 . CpI . Tx . dx 

A-0 

« a . D1 . CpI . (Til + Ton . 1/2 
Total heat content for n tubes is given by 

HP - n . a . D1 . CpI . 1 . (Til + Tol ) / 2 

or. 


HP ■ d 4 (Til + Tot ) 



£1 


Therefore rate of accumulat ion of heat in primary coolant tubes in term 
of deviation variables is given by 


d CHP3 d T i 1 d Tol 

_ - d4 . + aA ~ 3.8 

dt dt dt 


SECONDARY SIDE : Rate of heat accumulation is given by t he equation 

d EHS3 d 

[ A h Di£ HI + A (L - h) DoS C X H* + (1 - X) HI J ] 

dt dt 


Differentiating right hand side after substituting for X.Do2 from eqn. 
(A £ ) and t hen subst itut ing for X from equat ion ( A 3) in the expression 
we get 


d CHS3 d h 

« [ A Di£ HI - A Hs Do2S XS - A HI (1 - XS) Do£ 3 

dt dt 


d DoS 

+ l A (L - hS) Do£S k9 £ Hs / Di£ - HI / Ds33 

dt 


Using deviation variables above equation will be 


d CHS3 d Jh d DoS 

* bl . + b£ 3.8a 

dt dt dt 


Now combining equation (3.7) , 

d jh d Do2 t 

bf + b£ 

dt dt £ 

Combining equation (3.6) , (A 4) 


(A 4) , CA 3) and (3.8a) we get 

C b6 Do2 + b8 To£ + b9 Ft + blO Fi£ 

+ bl 1 Til 3 3.9 

, (A 3) end 





22 


d HI d Tol 

d4 + d 4 

dt dt 



t b6 Do2 + b8 To2 + b9 M + blO Fi2 

+ bl 1 Tit 3 


d Tol 

or, - A ( 4 , 2 ) Dog + A(4,4) Tol + B(4,1) £1 + B<4,2) Fig 

dt 

+ D( 4, 1) Til 3.10 


separating derivative terras from the equations (3.4) and (3.9) , that 
is , doing (3.4).b2 - (3.9).a2 we get 

d £ 

. A(1, 2) Do2 + A(1, 4) Tol^ + b(1,1) FI + B(1,2) Fi2 

dt 

+ DM,1) Ti_1^ 3.11 

and doing (3.9).a1 - (3.4) .bl we get 


d Do2 

« A(2, 2) Do2 + A(2, 4) Tol + B ( 2, 1 ) FI + B(2,2) Fig 

dt 

+ D(1 ,2) Til 3.12 


3.4-B STATE SPACE REPRESENTATION 

a) STATE SPACE EQUATIONS : Let the derivative ot a variable(say h) at 
tine k is 


h h 

d h k + 1 k 

dt k t 

where t is snail discretized tine 

Using the above def inat ion and equation (3.11) , (3.12) , (3.2) and 
(3.10) the state space raatrix eqation of the nodel is 
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h 


1 . 0 

A ( 1 / 2 ) . t 

0.0 

A(1 ,4) .t 


h 

DoS 

Af 


0.0 

0.0 

1 +A(£,£) , t 

0.0 

0.0 

0.0 

A (£, 4 ) . t 
0.0 


DoS 

Iol 

k + 1 

0.0 

A(4, £) . t 

0.0 

1 +A ( 4, 4 ) . t 

• 



B( 1 , 1 ) .t 

B ( 1 , £ ) .t 




D(1 ,1 ) .t 

+ 

B(£, 1 ) .t 

B ( 2 , £ ) . t 


FI 

| + 

D(£,1) .t 


c9 

0.0 


Fi2 


0.0 


B(4, 1 ) .t 

B(4,£) .t 



k 

D ( 4, 1 ) .t 

X « 

A X + 

B U + I 

D 

U 



- k + 1 

- - k 

- - k 

- k 
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c) MESUREHENT EQUATIONS 


y 


H X + V 

- - k 


3.14 


where y 

- k 


"iL 

AP 


k 


H 


r~ 

1 

0 


0 

0 


0 

1 


0 

0 


V * noise in the Measurements 
- k 


d) CONTROL LAU 


U 

- k 



3.15 


where C 


0.0 @1 

|S 0.0 


3.4-C STABILITY ANALYSIS 


Combining state space equation with the control law and measurement 


equation we get 



X 


■ (A +B. C . H) X + DU 
-k + t -k - - k 

If noise equals to zero , we have 


X » (A + B.C.H) X 

- 1 - 0 


£ 

X - (A + B.C.H) X * (A + B.C.H) X 

-£ - 1 - - 0 


k + 1 

X - (A + B.C.H) . X 

-K+1 - 0 


From the above equation it is obvious that the system will be stable 
iff (A + B.C.H) matrix have eigen values less then unity . Thus we 
have to find out suitable values of §1 , @£ and t , that will make 
the system stable ,ie, (A + B.C.H) matrix have eigen values less then 
unity . 

a) EXPRESSION FOR (A + B.C.H) MATRIX 

(A + B.C.H ) « 


1 +B 1 1 ,£).§£. t 
B t £, £ ) . §£ . t 
0.0 

Bt4,£) .|£.t 


At 1 ,£) .t 
1 +A (£, £ ) . t 
0 . 0 

A (4 , £ ) .t 


B( 1 , 1 ) . §1 . t 
B ( £, 1 ) .§1 . t 
c9 . §1 

B(4, 1 ) .|1 .t 


AC 1 ,4) .t 
At £, 4 ) - 1 
0.0 

1 +At 4/ 4 ) .t 


b) A computer programe is used to compute eigen values. Programs logic 
is shown in fig 3 . £ . 


c) RESULT OF STABILITY ANALYSIS : System will be in the most stable 
condition when maximum absolute value of eigen values is far less then 



unity . Fro* the computed result it is seen that the Host stable con 
dition is : 


t *1.0 sec 
ftl - 0 . 968E-08 
|£ - 0 . 239E-03 

3.4-D LEAK MODEL -: 

Loss of primary coolant through a leak or break within the primary cool- 
ing system (PHT-loop) is called loss of coolant. If liquid injection in- 
to the system is insufficient to cope up with this problem , it will 
cause an accident, called loss of coolant accident (LOCA). Depending up 
on the size of the leak it is classified into three categories ; 

- Small break LOCA ( leak die < 0.005 m) 

- Medium LOCA ( .005<leak dia < .025 m) 

- Large break LOCA ( leak dia > . 025 m ) 

The probabi 1 ity of occurrence of such a leak is very small but once it 
occurs it has high risk because it's consequence could include the risk 
of major release of radioactive materials . So PHT loop should be desi- 
gned aganist such fault . Thus it is a designed basis accident (DBA) . 
Table 3.T1 (ref. 6 Pp 292) shows the leak size of different LOCA, their 
probability of occurrence, percentage of fluid loss and no of injecting 
system should be operated to avoid bad consequence . 

From the table it is obvious that the occurrence of SBLOCA is 
higher and our consideration is restricted to such a leak . When such a 
leak occurs in a tube it has several consequence events depending upon 
size , tube size , orientation of the tube ( horizontal or vertical ), 
pressure difference at the escaping point , velocity of the fluid and 
nature of the fluid . Experimentally it is seen that phase separat ion 



will be Minimum when leak size is small , tube is vertical , pressure 
difference at the escaping point is small , velocity of fluid is high 
and fluid is highly viscous (ref. 3 Pp 15 - 19) . In our case , all the 
points except last one are in favour of no phase separat ion . Thus we can 
assume that no leak equations are equaly valid for leak model except 
the flow rate of primary coolant , FI , should be (FI - b) . where b is 
the magnitude of the leak . 

Now the state space representation of the LEAK MODEL is ; 

X « A . X + B.U + D.W + b . m 3.16 


k+1 - k - 

- k 

- k 


where the leak matrix. 

"bm.d.T 


~. 1 08E+00 - 

m * 

B(£, 1 ) .t 

as 

. 2S5E+0S 


0.0 

B( 4 . 1 ) .t 

i 

-. 131E+06 
- . £1 7E+02 


3.5 NUMERICAL VALUES OF THE MATRICES 


STATE TRANSITION MATRIX ( A ) 


. 001 E+0 
. OOOE+O 
. 000E+0 
.OOOE+O 


. 1S0E-3 
. 654E+0 
. 000E+0 
. 567E-S 


. OOOE+O 
.OOOE+O 
. OOOE+O 
. 000E+0 


-.838E-3 
- . 1 75E+0 
.000E+0 
. 831 E+0 


MEASUREMENT MATRIX ( H ) : 


. 1 00E+0 .000E+0 .000E+0 

.000E+0 .000E+0 .000E+0 


.OOOE+O 

.000E+0 


CONTROL MATRIX ( B ) PROCESS NOISE MATRIX ( D ) 


. 1 08E+00 . 331 E+0 


. 838E-3 

- . ££5E+2 .649E+2 


- . 1 75E+0 

. 1 31 E+6 .000E+0 


.000E+0 

. £1 7E+2 . 745E+£ 

1 

. 1 £4e+1 




READ INITIAL 
AND FINAL 
@1 & §2 


1 

COMPUTE (A+B.C.H) 

MATRIX AND EICEN 
VALUES 





WRITE §1 , §2 
AND EICEN 
VALUES 


FIG 3.2 : FLOW CHART OF STABILITY ANALYSIS 





TABLE - 3.T1 
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# CHARACTERISTICS OF LOCA AND FUNCTION OF YHE SAFETY SYSTEM « 



Large break 

LOCA 

Intermediate 

LOCA 

small break 

LOCA 

Break size (H 2 ) 

T 

Ul 

1 

OJ 

1 

Ul 

8E-3 - 4E-£ 

< 8E-3 

Ave . probability 




of occurrence 

IE-4 / a 

8E-4 / a 

3E-3 / a 

X loss of coolant 

10 - 15 

£ - 10 

< £ 

Safety system 

function (injection) 

- Accumulator 

3 of 4 

£ of 4 


- Auxiliary feed 

water 

— 



1 or £ of 4 




CHAPTER 4 


RESULTS AND DISCUSSION 


4.1 RESULTS 

The simulation has been done over the first 500 seconds 
of the accident , during which the operator has yet taken 
no corrective action . It is a real time simulation which 
predicts dynamic response to both the normal and abnormal 
operating conditions . The process model as given by the 
equations (3.13 - 3.16) are used for simulat ion . Since 
all the variables in the process model are deviation vari- 
ables , ie ( true value - normal value ) , henceforth in 
the discussion all variables mentioned are deviation vari- 
ables . 

At first the process is simulated without a leak to get 
the system dynamics . The procedure adopted is as follows: 

- initialize the state and control variables . 

- compute process noise by generating random 
error from multivariate normal distribution 
with zero mean and coverience matrix R 


- compute state variables by equation (3.13) 
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- generate measurement noise by sane procedure 
as in second step . 

- compute measurement variables by equation 
(3. 15) 

- use kalman filter to estimate the state vector 
and error coverience matrix . 

- obtain control inputs at each time using equa- 
t ion (3.14) 

The true values of the measured variables are obtained 
for 500 time steps , ie for 500 seconds , since each time 
step equals to one second . Fig. 4.1 shows the percent devi- 
ation of the measured variable pressure with t ime , which 
clearely indicates that the values are randomly distributed 
around zero mean ( ie normal value ) . 

Now for the same t ype of simulat ion a leak of different 
magnitude is incorporated in the dynamic model . Since the 
objective is to detect SBLOCA ( ie manitude of the leak is 
less than 2. OX of total flow rate ) , the leak magnitude is 
varied between 0.01 and 2 . OX . Simulation is done for 100 
times for each leak size . Table 4.T1 shows the magnitude of 
different leaks and accuracy with which it can be detected . 

From the table it can be seen that the method can de- 
tect leak of magnitude as low as 0.25X of total flow rate 
with a probability of 0.9 ( ie 90 times out of 100 times ) . 

For a leak of magnitude 0.1X , the percent deviation of pre- 
ssure is shown in the fig. 4.2 . From the figure it is clear 
that the percent deviation of total allowed deviation ( maxi- 
mum deviation which can be controlled ) is varying around 
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zero . This is the reason for not detecting this magnitude 
of the leak . Fig. 4.3 shows the percent deviation of pre- 
ssure for a leak of magnitude 1 . OX of total flow rate . The 
figure shows that the pressure will drop abruptly after 4.0 
seconds and then starts oscillating around 20. OX of the total 
allowed deviation . Thus the GLR method can detect this magni- 
tude of the leak . 

Fig. 4.4 shows the percent deviation of pressure for diffe- 
rent magnitude of leaks ; and fig. 4.5 shows the percent devi- 
ation of SG-water level . Both the figures show that the occu- 
rrence of leak effects the variables upto certain period of 
time and then the variables start oscillating around some 
"pseudo steady state" value . Thus the dynamics is similar to 
the response of a second order "single input-single output" 
(SISO) system to the step input . Now , if the time lapses by 
the system during which the variable should change by 66. 6% 
of total change is response time then for pressure drop it is 
5.0 seconds in all cases . Thus the response time for the pre- 
ssure drop is more or less constant . This can be explained as 
follows ; the occurrence of leak can be viewed as step input to 
the system . Now the response time will depend on the time that 
should be taken by the occurrence of leak to effect the flow 
rate of heavy water at the outlet of the U-tube . Since magni- 
tude of different leaks are small compared to total flow rate of 
heavy water , the velocity of water is more or less constant and 
it will take equal time to effect the outlet flow rate of the 
primary coolant ( heavy water ) at the outlet of the U-tube and 
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consequently the response time will be seme . 

Response time for pressure is 5.0 seconds which is very smell 
but compereble to the propegetion time from inlet to the outlet 
es explain below ; 

2 . 

totel cross scctionel eree of U-tubes(n . a )«1 950*8 . 92RE-5 M 
length of U-tubes ( 1 ) * IS. 8 H 

mess flow rete of heavy water ( ml ) ■ 380. SI Kg / S 
density of heavy water ( D 1 ) * 1100 Kg / 

Now time taken by heavy water to reach to the outlet of U-tube 
starting from inlet of the U-tube 

1 . n . a . 01 

SB 

ml 


IS. 8 * ( 1950 * 8.9SE-5 ) * 1100 
__ 


« 6.4 S 

If it is assumed that the occurrence of the leak will effect 
the system only when it's effect reaches to the outlet of the 
U-tube , expected response time will be 6.4 seconds . This is 
in close aggrement with the predicted value ( 5 . 0 seconds ) . 

From the figure ( fig. 4 . 4 ) it is obvious that the higher 
magnitude of leak , which will be viewed as a step input » 
effects the system to higher extent , ie the "pseudo steady 
state” will be far away from original steady state value . 
This is very similar to the response of a second order SISO 
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system to different step inputs . 

Fig. 4.5 shows that the response time for water level which 
depends upon the magnitude of the leak and in fact , higher the 
magnitude of the leak lower will be the response time . This is 
because higher magnutude of leak will effects the water level 
quickly due to radical change of rate of heat transfer and have 
low response time . Other things like effect of magnitude of leak 
on the new "pseudo steady state" can be explained as before . The 
response time for change of SG-water level is 25.0 seconds 
( Fig 4.5 ) . 

From the Fig. 4.4 A Fig. 4.5 one can notice that the res- 
ponse time for level ( 25.0 seconds ) is higher than the response 
time of pressure ( 5.0 seconds ) . This is because the occurrence 
of leak has direct effect on pessure and it has secondary ( indi- 
rect ) effect on the level . The response time are very low com- 
pare to the response time of most of the chemical plants . This 
is because of very high velocity of the primary coolant through 
the U-tubes . 


4.2 SCOPE OF THE MODEL AND THE METHOD 


Regarding the steam generator model , a number of questions 
can be considered more closely . Is the model accurate enough 
to represent the real situat ion or does it need a more accurate 
model and if so , what is the cost of computation ? Nonthless , 
the model as it is presented is quite simple , computationally 
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fast and gives satisfactory results in the steady state condi- 
tions and for the case of SBLOCA ( magnitude of the leak less 
than 2. OX of the total coolant flow rate ) . Thus such a model 
can be used with resonable confidence in numerous applications, 
such as steam generator primary tube break simulation . 

Regarding the performance of the method ( CLR ) , two impor- 
tant things should be considered simultaneously , the accuracy 
with which it can be detect SBLOCA and secondly , the delay 
time ie the time delay to detect the leak after the occurrence 
of the leak . The performance of the method is quite good be- 
cause it can detect a leak of magnitude as low as 0.25% of total 
flow rate with 90% accuracy . And the delay time is only 10-15 
seconds . 

Finally it can be concluded that the method can be used satis- 
factorily along with the suggested simplifications to detect 
SBLOCA . 

4.3 SCOPE OF THE WORK 


The model is not without its drawbacks . The primary disad- 
vantage to the model is an inability to handle two phase condi- 
tions within the primary loop . Attempts to remedy this short- 
coming should be made . Thus one can study the SC unit with a 
more realistic model . Another direction of the work is to simu- 
late the complete loop . For this , one needs to model the inter- 
mediate feed pump and pipe lines and couple them together to get 
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a consolidated model which can detect the leak in the entire PHT 
loop . 

In the design of the degital controller it is assumed that it 
follows linear control law , which is not always true , ie the 
transition from one steady state to another steady state . The 
CLR technique will not work for transition of the plant from one 
steady state to another . Effect of all states variables on the 
control law also can be verified . 

So far the CLR method is used in nuclear industry to detect 
the sensor failures and biases . Here in this work, it has been 
shown that the CLR technique can successfully be used to detect 
LOCA specially SBLOCA in the SC of the PHT loop . The detection 
of SBLOCA is critical as it may lead to an accident . Also 
mechanically its detection is difficult . The use of CLR method 
is illustrated in this respect and also got satisfactory perfor- 
mance of the method in detecting the SBLOCA . 

The method applied here can detect one fault at a time . 
Methods can can be devoloped to identify multiple faults (ref-18) 
of any type . Some other faults of greater importance in a nu- 
clear plant can be modelled such as coolant pump trip , abnormal 
valve beheviour , turbine trip auxiliary feed water failures , 
etc . Thus one can study the performance of the method in actu- 
al situat ions . 

This technique can also be developed for the actual measure- 
ment data instead of the innovations . For this further detailed 
study is needed >. Finally an expert system to detect SBLOCA can 
be developed combining the results of the CLR method with experi- 
mental knowledge as well as information from other sources . 
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TABLE -4.T1 


# # PERFORMA NCE OF THE CLR METHOD * * 


Magnitude of 

Leek in % 

NO. of tines Detected 

Out of 100 Tines 

% Accuracy 

, * 

0.01 

1 

0. 0£ 

0.10 

6 

0.06 

0 . 25 

90 

0.90 

0.50 

91 

0.91 

1 . 00 

96 

0.96 

1 .50 

98 

0.98 

£.00 

98 

0.98 





OF TOTAL ALLOWED 
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APPENDIX 


A . 1 CORRELATION BETEUEEN X 6 Do£ 

From the defination of steam quality we have ; 
•ass flow rate of dry steam * steam quality X 


ie, 

mg - 

X 

. FoE 

. DoS 

and 

ml ■ 

(1 

- X) . 

FoS. Do£ 


Total volumetric flow rate * flow rate of steam 


i 

Fo£ « mg / Ds 

+ ml / DiS 

combining above 

three 

equat ions we get 


1 / Do£ 

- X / 

Ds + (1 - X) / DiS 

or. 

V 

Ds . 

(DiS - DoS) 

A m 

DoS. 

(DiE - Ds) 

and. 

X.DoE « 

( 

DiE - DoS) . Ds 


( DiE - Ds) 


Writing (A 1 ) in term of deviation variables we get 


( X + XS ) 


Di£ . Ds 

+ 

(DIE - Ds) . (Do£S + Do£ ) 


Di£ . Ds 


(Di£ - Ds ) . Do£S 


Do£ 

1 + 

Does 


DiS . Ds 
(Di£ - Ds) .DoSS 


DoS 


1 

DoES 


(total mass flow rate 
of vapour-water mixture 


+ water 


(At ) 


(AS) 


Ds 


(DiS - Ds) 


Ds 


(DiE - Ds) 


Ds 


(Die - Ds) 


( neglecting higher terms ) 



At steady state deviations are zero , thus from the above equation 
we get 


Die . Ds Ds 

(Di2 - Ds) .Does (Die - Ds) 

Combining above two equations we get 


Die . Os 

X - - . Doe 

(Die - DsO) .Does 

or, X - k9 . Doe (A3) 

A.e EXPRESSION FOR THERMAL IMBALANCE IN THE STEAM GENERATOR 

Thermal imbalance at any time in the SG is defined as the difference 
of heat content between incoming and outgoing streames ie 

B ■ Rate of total heat in - rate of total heat out 

« (FiE. DiS. HI + FI .01 .CpI .Til ) - (Foe . DoS . CX. Hs + (I-X)HIJ + 

FI . D1 .CpI .To! ) 

Using deviation variables and equqt ions (A 1) & (A e> we get 

_B - -Foes. Does. (Hs - Hl).k9. DoS - FIS. Dl. CpI. Tot -D1 . CpI . ( Ti 1 S 
- TotS) .£1 - Foes (Hs - HD.XS.DoS + DiS. HI. £ie + F1S.D1.Cpt. Til 

or, B - b6 . DoS + b8. Tot + b9. £1 + blO. FiB + bll. HI (A4) 

A. 3 CALCULATION OF THE ELEMENTS OF MATRICES 

Nuclear power plant can be operated at any one of the following 
stages ; 


- steady state operation with full load 

- steady state operation with reduced load 



Our aodcl can be used for both the case* . But here only the former case 


is analyse , ie , the elements of matrices are calculated for steady 
state condition . Table A- 1 shows the value of variables at steady state 
condition . And table A-£ contains values of constant parameters . Diffe- 
rent constant coefficients of the modelled equations are calculated based 
on these tables . Calculation is as follows . 

a) / 

al - (A . Di£ - A . Does) 

- (1.644 * 852.5 - 3.973 * 18.60e> 

- 1387.79 Kg / n 
a2 » A. (L - hS) 

- 3.973 * (8.0 - 6.4) 

« 6.356 H 3 

a6 - - Fo2S * -8.04 M / S 
al 0 * Di2 

« 858.5 Kg / M 3 
/ 

bl - A . Di8 . HI - A . Do8S CHs . XS + (1 - XS) . HI) 

- 1.644 * 852.5 * 1085.36 - 3.973 * 18.608 C 28o1 .5 * .9 

+ Cl -.9) * 1085.363 

« 1386996.08 KJ j H 

b8 * A . ( L - hS) . DoSS . k9 . ( Hs / Di8 - HI / Ds ) 

■ 3.973 * ( 8.0 - 6.4 ) * 18.602 * (- 0.04946) * 

( 2801.5 / 858.5 - 1085.36 / 16.779) 

- 359.08 KJ . M 3 / Kg 

b6 « - ( Fo2S . DoSS ) . ( Hs - HI ) . k9 - FoSS.XSCHs - HI) 

■ - 38.021 * ( 2801.5 - 1085.36) * ( - 0.04946 ) 

3 - 2. 04*. 9(2801 .5 - 1085.4) 

-70.41 KJ . My Kg . S 

b8 - - ( FIS . D1 ) . CpI 

- - 380.21 * 5.518 

« - 2097.84 KJ / S . k 



b9 * D1 . CpI . ( TilS - TolS) 


- 5.518 * 1100 * ( £93.4 - 249) 

- £69479.6 KJ j M 3 

blO - Di£ . HI 

« 852.5 * 1085.36 
■ 925269.4 KJ j M 3 

bll - ( FIS . D1 ). CpI 
« 380.21 8 * 5.518 
« £097.84 KJ . S j k 

d4 ■ n . 1 . a . CpI / 2 

- 1950 * 12.8 * 8 . 1 97E-05 * 5.518 ♦ 1100 

- 6208.82 KJ / k 

k 9 •= - Di2 . Ds / { ( Di£ - Ds ) . Do2S ) 

- - 852.5 * 16.779 / C ( 852.5 - 16.779 ) * 18.602 

« - 0.04946 If 3 / Kg 

k 1 0 - *1 . b2 - bl . *2 

■ 1327.79 * 359.08 - 1326996.02 * 6.357 

- - 7957737.5 KJ / H 

A( 1 ,2) - - b6 . *2 / ( 2. 0 . ktO ) 

- - 70.41 * 6.356 / C £.0 * ( - 7957737.5 ) J 

■ 1.201 5E-04 mV Kg . S 

A( 1 ,4) * - b8 . «2 / ( 2.0 . k 1 0 ) 

« - 2097.8 * 6.356 / C 2.0 * ( - 7957737.5 ) 3 

- - 8 . 378E-04 M / S . k 



A-5 


A(£, 2) - b6 . *1 / ( £.0 . klO 0 

- 70.41 * 1327.8 / C £.0 ( - 79957737.5 ) 3 

- - 0.346 / S 

A( 2, 4 ) - b8 . *1 / ( £.0 . k 1 0 ) 

« £097.8 * 1327.8 / C £.0 ( - 79957737.5 ) 3 
■ - 0.175 M/S k 

A (4,2) - b6 / ( £.0 . d4 ) 

- 70.41 / ( £.0 * 6208.8 ) 

- 0 . 567E-02 M 3 / Kg . S 

A( 4 , 4 ) « b8 / ( 2.0 . d4 ) 

- - £097.8 / ( £.0 * 6£08 . 8 ) 

- -0.169 /S 

B( 1 , 1 ) ■ - b9 . *2 / ( £.0 . k 1 0 ) 

- - 869479.6 * 6 . 356 / C £ . 0 ( - 7957737.5 ) 3 

- 0.108 /H Z 

B( 1 ,£) * ( *10 . b£ - blO . *2 / £.0 ) / klO 

- ( 852.5 * 359.08 - 925269.4 + 6.356 / £.0 ) / ( - 7957737.5) 

- 0.331 /M 2 

B(£, 1 ) - b9 . *1 / ( 2- 0 . k 1 0 ) 

- 269479.6 * 1327.8 / C 2.0 ( - 7957737.5 ) 3 

✓ 6 

- - 22.48 Kg / M 

B(2,2) * ( blO . *1 / £.0 - *10 . bl ) / klO 

■ (925269.4 * 1327.8 / £ - 852.5 * 1326996.01 ) / (-7957737.5) 

- 64.883 .Kg / M* 

B(4, 1 ) - b9 / ( 2. 0 . d4 ) 

- 269479.6 / ( 2.0 * 6208.8 ) 

• 21.7 k / M 3 



A-6 


8(4, £) - blO / ( £.0 . d4 ) 

- 925269.4 / ( £.0 * ££08.6 ) 

- 74 . 51 £ k / M 3 

0(1, 1 ) - - bll . a£ / ( £.0 . k 1 0 ) 

- - £097.8 * 6.356 / C 2.0 ( - 7957737.5 ) 

- 8 . 378E-04 M / S k 

D(£, 1 ) « bll .al / ( £.0 . k 1 0 ) 

- £097.8 * 13£7.8 / £2.0 ( - 7957737.5 ) 1 

- -0. 175 Kg / M 3 . 8 . k 

0(4, 1 ) « bl 1 / d4 - 1 .414 / t 

- £097.8 / ( £.0 * 6£08 . 8) - 1.414 / t 

- ( 0.169 - 1.414 / t ) « - 1 .£45 /S ( sine* t-1 . 0) 


b) CALCULATION OF c9 ( r*f. 13 Pp £54 ) 


Reynold's no N 

re 


V . id . 01 id . ( FIS . 01 ) 


( mu ) ( mu ) . £ £ 

n a 


0.0102 * 380.21 

m mm, mm mm mm mm mm mm mm mm mm mm mm mm mm mm mm mm mm mm 

10E-0S * 1950 * 8 . 1 97E-05 

= 2.43E+06 


At this N friction factor , f - 0.007 ( ref. ££ Pp £54 > 

re 


thus using this value of f we get , 

c9 « f . 1 . ( FIS . 01 ) / ( id . n .a ) 

- 0.007 * 1 £ . 8 * 380.21 / ( .0102 * 1950 * 8.197 ) 

4 

■ 130519.8 Kg / S . 



TABLE - A.T1 
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«tt LIST OF THE VARIABLES AND THEIR STEADY STATE VALUES ## 

1 . hS - height of the liquid level in SC : 6.4 M 


2. Do2S - density of outgoing steam 

3. PS - pressure drop in primary coolant tubes 

4. TolS - outlet temperature of primary coolant 

5. F1S.D1 - mass flow rate of primary coolant 

6. Fo2S.Do£S - mass flow rate of secondary water 

7. XS - steam quality 


16.602 

Kg / M 

1 3E+0E 

£ 

N / M 

522 k 


380.21 

Kg / S 

38.021 

Kg / S 

0.9 



TABLE - A.T2 


## LIST OF THE PARAMETERS AND THEIR VALUES #ft 

3 


1 . 

D1 

— 

average density of water 

1 

1100 Kg / M 

2. 

Ds 

- 

density of saturated steam at 521 k 

: 

3 

16.779 Kg / M 

3. 

Di2 


density of saturated water at 521 k 


3 

852.5 Kg / M 

4 . 

Til 

- 

inlet temperature of primary coolant 


566.4 k 

5. 

TB 

- 

temperature of the feed water 


521 k 

6. 

HI 

- 

enthalpy of saturated water at 521 k 


1085.36 KJ / Kg 

7. 

Hs 

- 

enthalpy of saturated steam at 521 k 


2801 .5 KJ / Kg 

8. 

CpI 

- 

average specific heat of heavy water 


5.518 KJ/Kg/ k 

9. 

mu 


viscosity of heavy water at 523 k 


1 . OE-5 Kg/M/S 

10. 

a 

- 

cross sectinal area of a coolant tube 

• 

£ 

8.197E-S M 

1 1 . 

A 

- 

cross sectional area of the upper 

section of the steam generator 

: 

3 . 973 M 2 

12. 

A 


cross sectional area of the lower section 
of the steam generator (excluding tubes) 

: 

1.644 M 

13. 

n 

- 

no of tubes per one steam grnerator 

: 

1950 

14. 

1 

- 

length of the U-tubes 

: 

12.8 M 

15. 

id 

- 

internal diameter of a single U-tube 

: 

1.02E-2 M 

16. 

L 

- 

shell height 

: 

8.0 M 



